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1.  Summary 

For  the  past  three  years,  we  have  been  developing  an  Exact  Scientific  Computational  Library  (ESCL) 
using  p-adic  arithmetic.  New  algorithms  have  been  designed  and  implemented  for  matrix  operations  with 
rational  numbers  by  representing  numerator  and  denominator  with  arbitrary  length  integers,  all  integers 
and  fractional  numbers  are  represented  by  p-adic  sequences,  and  all  arithmetic  calculations  are  carried 
out  in  p-adic  domain.  In  this  project,  we  have  worked  on:  1)  investigating  the  relation  of  the  length  M  of 
p-adic  expansion  for  a  rational  matrix  and  the  periodicity  of  a  resulted  p-adic  sequence  from  arithmetic 
operation  in  p-adic  field;  and  extension  of  the  ESCL  to  compute:  2)  the  complex  rational  matrix,  3)  the 
exponential  of  a  rational  matrix. 


2.  Progress  on  Length  M  and  Periodicity  of  a  p-adic  Expansion 


To  determine  what  is  the  efficient  length  M  of  the  p-adic  expansion  for  a  rational  number  in  the  matrix 
operation  is  not  a  trivial  problem.  Let  us  observe  what  happens  after  the  arithmetic  operations  of  two  p-adic 
sequences. 

All  rational  numbers  can  be  uniquely  written  in  the  form, 

a  =  YjajPJ  ■  (D 
o 

We  know  that  a  real  number  is  rational  if  and  only  if  its  decimal  expansion  is  periodic.  Similarly,  a 
p-adic  number  is  rational  if  and  only  if  its  p-adic  expansion  is  periodic.  Consequently,  since  we  are 
primarily  interested  in  the  p-adic  expansions  of  rational  numbers,  we  will  be  dealing  only  with  p-adic 
expansions  which  are  periodic.  The  expansion  eventually  repeats  to  the  right.  That  is,  if  a  is  a  rational 

number,  then  it  has  a  repeating  pattern  of  a tpJ  in  its  p-adic  expansion,  i.e.,  it  is  of  the  form 

a  =  .A0...Asa0...a„_, .  (2) 


Addition/subtraction 

Assume  that  we  have  two  p-adic  sequences,  (s<t): 
a  =  .A,...Asal...an 
b  =  .Bt...B,b,...b„ 

Considering  various  carry  digits’  effects,  we  concluded  that  the  maximum  length  of  the  p-adic 
expansion  of  (a  +  b)  is: 

2xLCM(n,m)  +  max(s,t)-\ .  (3) 
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Multiplication 

axb  =  .Ai...Asal...a„x.B]...Blb]...bm 


=  .Ax...Asax...anx.Bx...B,  +  .Ax...Axax...anx.0...0bx...blf 


=  .Ax...Asax...a„x.Bx...B,+.Ax...  As  x  .0...0  bx  ...bm+  .0...0a,  ...anx  .0...0  bx .  ..bm 

v - - - v — - - J  v - - - v - - - J  v - v - ^ 

1  2  3 

Considering  all  the  three  parts  separately  and  the  various  carry  digits’  effects,  we  concluded  that,  in 
multiplication,  the  length  of  periodic  part  of  the  product  is: 


LCM(m,  n)  x  (pGCD^n’>  - 1) ,  (4) 

where  m  and  n  are  the  length  of  periodic  part  of  the  two  multipliers,  p  is  the  prime. 


The  length  of  periodicity  of  the  resulting  p-adic  sequence  can  be  very  large  from  (4).  But  if  we  should 
represent  all  the  p- adic  sequences  with  a  complete  period  during  all  the  calculations,  we  can  definitely  carry 
out  all  the  arithmetic  operations  exactly.  Further  investigation  is  needed  to  find  proper  length  M  which  is 
smaller  than  the  periodicity  of  a  p-adic  expansion.  The  determination  of  the  efficient  length  M  of  the  p-adic 
expansion  for  a  rational  number  in  the  matrix  operation  becomes  more  important,  when  the  algorithms  in 
eigensystem  computation  use  iterative  arithmetic  operations. 


3.  Progress  on  the  Computation  of  Complex  Rational  Matrices 

We  have  implemented  the  algorithms  developed  for  single  and  matrix  rational  number  operations 
(addition,  subtraction,  multiplication  and  division)  in  the  p-adic  field,  to  compute  the  complex  rational 
arithmetic  operations  in  the  p-adic  field  as  follows: 


Complex  addition 

(a  +  bi)  +  (c  +  di )  =  (a  +  c)  +  i(b  +  d) ,  (5) 

Complex  subtraction 

(( a  +  bi)  -  (c  +  di)  =  (a  -  c)  +  i(b  -  d) ,  (6) 

Complex  multiplication 

(a  +  bi){c  +  di)  =  ( ac  -  bd)  +  i(ad  +  be) ,  (7) 

and  Complex  division 

a  +  bi  ac  +  bd  .be-  ad 

- =  -; - T  +  l~i - 77-  (g) 

c  +  di  c2+d '  c  +d 


These  complex  operations  are  needed  in  the  process  of  eigenvalue  computation. 
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4.  Progress  on  the  Computation  of  Matrix  Exponential 

4.1  Compute  matrix  exponential  with  Jordan  canonical  form 

Matrix  decomposition  methods,  which  are  to  be  most  efficient  for  problems  involving  large  matrices 

evaluation  ofe'A  ,  are  those  which  are  based  on  factorizations  or  decompositions  of  the  matrix  A.  The 
Jordan  canonical  form  ( JCF)  decomposition  is  one  of  the  many  stated  in  [1].  Due  to  its  truncation  error  by 
floating-point  arithmetic  operations,  it  has  not  been  widely  used,  since  a  single  rounding  error  may  cause 
multiple  eigenvalues  to  become  distinct,  which  can  alter  the  entire  structure  of  the  decomposition.  If  we  can 
improve  the  accumulated  floating-point  truncation  error  created  during  the  iteration  process  of  the 
computation  by  p-adic  arithmetic  wherever  possible,  we  hope  that  the  stability  and  accuracy  can  be 
improved  dramatically. 


Definition:  The  Jordan  canonical  form  decomposition  states  that  there  exists  an  invertible  P  such  that 
P~x  AP  =  J ,  (9) 

where  J  is  a  direct  sum,  J  =  J,  0  •  •  •  ®  Jk ,  of  Jordan  blocks, 


A,  1  0  •••  0 

0  A,  1  •••  0 

1 

0  0  0  •••  A, 


(m,  -by-m,). 


(10) 


The  At  are  eigenvalues  of  A,  here  we  assume  that  A  is  rational.  Each  Jordan  block  corresponds  to  a 


linearly  independent  eigenvector.  If  any  of  the  mj  is  greater  than  /,  A  is  said  to  be  defective.  This  means 
that  A  does  not  have  a  full  set  of  n  linearly  independent  eigenvectors.  A  is  derogatory  if  there  is  more  than 
one  Jordan  block  associated  with  a  given  eigenvalue. 

For  a  given  matrix  A,  its  Jordan  canonical  form  J  is  completely  determined  by  the  maximal  number 
of  linearly  independent  eigenvectors  of  A:  the  number  of  the  Jordan  blocks  in  J  is  equal  to  the  maximal 
number  s  of  linearly  independent  eigenvectors  of  A,  each  of  which  is  associated  with  a  Jordan  block  whose 
order  is  the  same  as  the  ‘rank’  of  the  eigenvector. 

Thus  the  number  of  Jordan  blocks  with  the  same  eigenvalue  A  is  equal  to  the  dimension  of  the 
eigen-space  Ex  =  N(AI  -  A) ,  the  number  of  linearly  independent  eigenvectors  of  A  belonging  to  A  . 
Moreover,  the  sum  of  the  orders  of  all  Jordan  blocks  associated  with  an  eigenvalue  A  is  equal  to  the 
multiplicity  mx  of  A  . 

In  principle,  the  problem  posed  by  defective  eigensystems  can  be  solved  by  resorting  to  the  Jordan 
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canonical  form  (JCF). 
If 


A  =  P[JX®  ■@JIC]P-'  (11) 

is  the  JCF  of  A ,  then 

e'A  =  p[eu>  ®---®eUt  Jp-1 .  (12) 


Thus,  it  is  enough  to  compute  eu  for  a  simple  Jordan  block  J.  The  exponentials  of  the  Jordan  blocks 
Ji  can  be  given  in  closed  form.  For  example,  if 

U.  1  0 


J.  = 


then, 


1 

A. 


u. 

e  =e 


(13) 


1  t  t2/  2! 

0  1  / 

1 

0 


r2/(n-  2)! 


(14) 


It  is  important  to  know  how  sensitive  a  quantity  is  before  its  computation  is  attempted. 

6 


One  of  the  most  important  issues  in  the  computation  of  matrix  exponential  with  JCF  is  the  evolution 
of  performing  eigenvalue  analysis. 


4.2  Power  method 

Numerical  analysis,  at  its  simplest,  is  an  iterative  technique.  It  is  a  fruitful  exercise  to  study  the 
Power  method  to  get  an  in-depth  understanding  of  the  numerical  solution  for  eigen-values  and 
eigenvectors. 


4.2.1  Algorithm 

The  Power  method  is  an  iterative  technique  used  to  determine  the  dominant  eigenvalue  of  a 
matrix — that  is,  the  eigenvalue  with  the  largest  magnitude.  By  modifying  the  method  slightly,  it  can  also  be 
used  to  determine  other  eigen-values.  One  useful  feature  of  the  Power  method  is  that  it  produces  not  only 
an  eigenvalue,  but  also  an  associated  eigenvector.  In  fact,  the  Power  method  is  often  applies  to  find  an 
eigenvector  for  an  eigenvalue  that  is  determined  by  some  other  means. 

To  apply  the  Power  method,  we  assume  that  the  nxn  matrix  A  has  n  eigen-values  A],X2,---,An 
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with  an  associated  collection  of  linearly  independent  eigenvectors  |v(l),  V(2),-  v(/,)  |.  Moreover,  we 

assume  that  A  has  precisely  one  eigenvalue,  \  ,  that  is  largest  in  magnitude,  so  that 


|A1|>|A2|>|^|>->|An|>0. 

If x  is  any  vector  in  9?rt,  the  fact  that  |v(l),v<2), 

n 

constants  •••,/?„  exist  with  X  =  ^/?yV0). 


>=i 


Mim) 


By  making  a  series  of  calculations,  we  get 


( j ) 
Pm- 1 


/K’, 


(j) 
Pm- 1 


,  v(w) j  is  linearly  independent  implies  that 


(15) 


and 


*(m)  = 


™  „(0) 


(16) 


By  examining  Eq.  (15),  we  see  that  limm_>00  p(m)  =  A, .  Moreover,  the  sequence  of  vectors  {x(m)  }“=0 


converges  to  an  eigenvector  associated  with  A, . 

The  Power  method  implementation  can  be  stated  as  follows: 

INPUT:  dimension  «;  matrix  A\  vector  x;  tolerance  TOL;  maximum  number  of  iterations  N. 


OUTPUT:  approximate  eigenvalue  p  ;  approximate  eigenvector  x  (with  ||x||  =  1 )  or  a  message  that  the 

maximum  number  of  iterations  was  exceeded. 

Step  1  Setk=l. 

Step  2  Find  the  smallest  integer  P  with  1  <p<n  and  )■*,,  |  =  • 

Step  3  Set  x-x/xp. 

Step  4  While  (k  <  N)  do  Step  5-11. 

Step  5  Set  y  =  Ax  . 

Step  6  Set  p  =  yp. 

Step  7  Find  the  smallest  integer  p  with  1  <  p  <  n  and  |>>r|  =  ||.p||oo . 

Step  8  If  yp  =  0  then  OUTPUT  (‘Eigenvector’,  x); 

OUTPUT  (‘A  has  the  eigenvalue  0,  select  a  new  vector  x  and  restart’); 


STOP. 


step  9  Set  ERR  =  \x-(y/yp)jj 

x=y/yP 

Step  10  If  ERR  <  TOL  then  OUTPUT  ( /u,x ); 
(The  procedure  was  successful.) 

STOP. 

Step  1 1  Set  k=k+ 1 . 


Step  12  OUTPUT  (‘The  maximum  number  of  iterations  exceeded’); 
(The  procedure  was  unsuccessful.) 


STOP. 


As  we  can  see,  such  an  algorithm  is  designed  to  get  the  approximate  eigen-values,  round-off  and 
truncation  errors  may  be  introduced  by  using  floating  point  arithmetic  during  the  computation,  but  the  exact 
linear  computation  can  avoid  that  kind  of  errors. 

4.2.2  Implementation 

For  example,  the  matrix 


-4 

14 

-5 

13 

-1 

0 

0 

0 

2 


has  eigen-values  \  =  6  ,  A2  =  3 ,  and  =  2  ,  so  the  Power  method  described  in  this  algorithm  will 


converge.  Let  x(0)  =  (1,1, 1)* ,  then 


y(,)=Ax(0)  =(10,8,1)', 


so 

vo> 

l/i  =10,//(,)=^0)=10,  and  X(1)  =  =  (1,°.8,0.1)' . 

II  lloo  10 

Continuing  in  this  manner  leads  to  the  approximations  to  the  dominant  eigenvalue  1/3. 

Comparison  of  the  results  of  the  exact  linear  computation  and  results  using  floating  point  arithmetic. 
The  results  using  floating  point  arithmetic  as  follows: 


1 

0.55555555555555558 

8 

0.3340807174887892 

2 

0.40000000000000002 

9 

0.33370618941088748 

3 

0.36111111111111094 

10 

0.3335195530726256 

4 

0.34615384615384615 

11 

0.3334263912153359 

5 

0.33950617283950624 

12 

0.33337984928830566 

6 

0.3363636363636362 

13 

0.33335658806567126 

7 

0.33483483483483478 

14 

0.33334495988838486 
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15  , 

0.33333914640810103 

38 

0.33333333333402604 

16 

0.33333623982003013 

39 

0.33333333333367965 

17 

0.33333478656401028 

40 

0.33333333333350645 

18 

0.3333340599455038 

41 

0.33333333333342008 

19 

0.33333369663862655 

42 

0.33333333333337667 

20 

0.33333351498578179 

43 

0.33333333333335502 

21 

0.33333342415950795 

44 

0.33333333333334414 

22 

0.33333337874640823 

45 

0.3333333333333387 

23 

0.33333335603986758 

46 

0.33333333333333592 

24 

0.33333334468659981 

47 

0.33333333333333437 

25 

0.33333333900996609 

48 

0.33333333333333359 

26 

0.33333333617164973 

49 

0.33333333333333326 

27 

0.33333333475249138 

50 

0.33333333333333326 

28 

0.33333333404291232 

51 

0.33333333333333326 

29 

0.33333333368812279 

52 

0.33333333333333326 

30 

0.33333333351072825 

53 

0.33333333333333326 

31 

0.33333333342203075 

54 

0.33333333333333326 

32 

0.33333333337768223 

55 

0.33333333333333326 

33 

0.33333333335550774 

56 

0.33333333333333326 

34 

0.33333333334442072 

57 

0.33333333333333326 

35 

0.33333333333887694 

58 

0.33333333333333326 

36 

0.33333333333610504 

59 

0.33333333333333326 

37 

0.33333333333471904 

60 

0.33333333333333326 

The  results  of  the  exact  linear  computation  using  p-adic  arithmetic  as  follows: 


1: 

7: 

5/9 

223/666 

2: 

8: 

2/5 

149/446 

3: 

9: 

13/36 

895/2682 

4: 

10: 

9/26 

597/1790 

5: 

11: 

55/162 

3583/10746 

6: 

12: 

37/110 

2389/7166 
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13: 

3758096383/11274289146 

14335/43002 

32: 

14: 

2505397589/7516192766 

9557/28670 

33: 

15: 

15032385535/45097156602 

57343/172026 

34: 

16: 

10021590357/30064771070 

38229/114686 

35: 

17: 

60129542143/180388626426 

229375/688122 

36: 

18: 

4008636 1 429/ 1 20259084286 

152917/458750 

37: 

19: 

2405 1 8 1 68575/72 1 554505722 

917503/2752506 

38: 

20: 

1 60345445717/481036337 1 50 

611669/1835006 

39: 

21: 

962072674303/28862 1 8022906 

3670015/11010042 

40: 

22: 

641381 782869/1924 1 45348606 

2446677/7340030 

41: 

23: 

38482906972 15/11 544872091642 

14680063/44040186 

42: 

24: 

2565527 1 3 1 477/769658 1 394430 

9786709/29360126 

43: 

25: 

153931 62788863/46 1 79488366586 

58720255/176160762 

44: 

26: 

1 0262 1 08525909/30786325577726 

39146837/117440510 

45: 

27: 

6 1 57265 1 1 55455/ 1 847 1 7953466362 

234881023/704643066 

46: 

28: 

4 1 048434 1 03637/123 1 453023 10910 

156587349/469762046 

47: 

29: 

24629060462 1 823/73887 1 8 1 3865466 

939524095/2818572282 

48: 

30: 

1 64 1 937364 1 4549/49258 1 209243646 

626349397/1879048190 

49: 

31: 

985 1 624 1 8487295/295548725546 1 882 
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50: 

656774945658 1 97/1970324836974590 

51: 

3940649673949 183/11821 94902 1 847546 

52: 

2627099782632789/7881299347898366 

53: 

15762598695796735/47287796087390202 

54: 

1 0508399 1305311 57/3 1 525 197391 593470 
55: 

63050394783 1 86943/ 1891511 84349560826 
56: 

42033596522124629/ 1 26 1 00789566373886 
57: 

25220 1 579 1 32747775/756604737398243322 
58: 

1 68 1 343860884985 1 7/504403 1 58265495550 
59: 

1 0088063 1 653099 1 1 03/30264 1 8949592973306 
60: 

672537544353994069/20 1 76 1 263306 1 982206 
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From  the  above  results,  we  can  find  that  floating-point  operation  does  not  function  very  well  in 
terms  of  round-off/truncation  errors.  The  approximation  of  the  dominant  eigenvalue  remain 
0.33333333333333326  from  loop  49,  while  it  should  be  0.33333333333333333  precisely.  The  results  have 
been  rounded  up  or  down  in  every  step,  finally,  the  iterative  process  will  dramatically  enhance  the 
uncertainty  of  the  outcome.  In  the  meantime,  using  p-adic  arithmetic  will  get  a  decisive  advantage  since  it 
can  keep  the  calculation  free  of  round-off/truncation  errors. 

An  obvious  advantage  of  the  exact  linear  computation  using  p-adic  arithmetic  is  that,  we  can  get 
results  as  accurate  as  we  want  with  the  loop  increasing,  while  floating-point  arithmetic  can  only  reach  a 
certain  precision  and  stay  at  a  value  no  matter  how  many  loops  it  runs. 

Compare  to  the  program  using  floating-point  arithmetic  in  which  the  results  remain 
0.33333333333333326  from  loop  49,  the  program  using p-adic  arithmetic  can  get  more  and  more  accurate 
after  each  loop,  for  example,  in  loop  59,  the  approximate  eigenvalue  is 

1008806316530991 103/3026418949592973306  *  0.33333333333333333366375685 , 

after  running  one  more  loop,  the  approximate  eigenvalue  is 

672537544353994069/2017612633061982206*0.33333333333333333349854509 

which  is  closer  to  1/3. 

Precision  is  critical  to  the  computation  of  the  exponential  of  a  rational  matrix  since  the  method 
requires  a  tremendous  amount  of  operations,  the  truncation  errors  for  every  step  will  cumulate  and  finally 
make  a  great  difference  from  the  exact  value.  As  a  matter  of  fact,  a  little  error  will  alter  the  entire  structure 
of  j  and  P  and  lead  to  an  utterly  wrong  outcome,  this  will  put  the  matter  beyond  a  doubt  that  our  “Exactly 
Computing”  system  would  prove  highly  valuable  on  the  computation  of  the  exponential  of  a  rational 
matrix. 

4.3  Eigenvalues  of  a  Real  General  Matrix 

The  solution  of  eigensystems  is  a  fairly  complicated  business,  almost  all  routines  in  use  nowadays 
trace  their  ancestry  back  to  routines  published  in  Wilkinson  and  Reinsch’s  Handbook  for  Automatic 
Computation,  Vol.  II,  Linear  Algebra  [2].  A  public-domain  implementation  of  the  Handbook  routines  in 
FORTRAN  is  the  EISPACK  set  of  programs  [3].  It  includes  the  ability  to  solve  for  eigenvalues  and 
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eigenvectors  of  various  kinds  of  matrices,  has  been  implemented.  The  routines  we  used  are  translations  of 
either  the  Handbook  or  EISPACK  routines. 

The  matrix  is  first  balanced.  Orthogonal  similarity  transformations  are  used  to  reduce  the  balanced 
matrix  to  a  real  upper  Hessenberg  matrix.  The  implicit  double-shifted  QR  algorithm  is  used  to  compute  the 
eigenvalues  and  eigenvectors  of  this  Hessenberg  matrix. 

4.3.1  Balancing 

The  idea  of  balancing  is  to  use  similarity  transformations  to  make  corresponding  rows  and  columns 
of  the  matrix  have  comparable  norms,  thus  reducing  the  overall  norm  of  the  matrix  while  leaving  the 
eigenvalues  unchanged.  It  is  recommended  to  always  balance  non  symmetric  matrices.  It  never  hurts,  and  it 
can  substantially  improve  the  accuracy  of  the  eigenvalues  computed  for  a  badly  balanced  matrix. 

Balancing  is  a  procedure  with  of  order  N2  operations.  The  actual  algorithm  used  is  due  to  Osborne, 
as  discussed  in  [2].  It  consists  of  a  sequence  of  similarity  transformations  by  diagonal  matrices.  The  output 
is  a  matrix  that  is  balanced  in  the  norm  given  by  summing  the  absolute  magnitudes  of  the  matrix  elements. 

Note  that  if  the  off-diagonal  elements  of  any  row  or  column  of  a  matrix  are  all  zero,  then  the  diagonal 
element  is  an  eigenvalue.  If  the  eigenvalue  happens  to  be  ill-conditioned  (sensitive  to  small  changes  in  the 
matrix  elements),  it  will  have  relatively  large  errors  when  determined  by  the  routine  hqr.  We  could  have 
determined  the  isolated  eigenvalue  exactly  and  then  deleted  the  corresponding  row  and  column  from  the 
matrix. 

4.3.2  Reduction  to  Hessenberg  Form 

First  we  reduce  the  matrix  to  a  simpler  form,  and  then  we  perform  an  iterative  procedure  on  the 
simplified  matrix.  The  simpler  structure  we  use  here  is  called  Hessenberg  form.  An  upper  Hessenberg 
matrix  has  zeros  everywhere  below  the  diagonal  except  for  the  first  sub-diagonal  row.  For  example,  in  the 
6*6  case,  the  nonzero  elements  are: 
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X  X  X  X  X  X 

X  X  X  X  X  X 

X  X  X  X  X 

X  X  X  X 

XXX 
X  X 

Here  we  use  a  procedure  analogous  to  Gaussian  elimination  with  pivoting.  Before  the 

rth  stage,  the  original  matrix  A  =  Ax  has  become  Ar ,  which  is  upper  Hessenberg  in  its  first  r-1  rows 

and  columns.  The  rth  stage  then  consists  of  the  following  sequence  of  operations: 

Find  the  element  of  maximum  magnitude  in  the  rth  column  below  the  diagonal.  If  it  is  zero,  skip 
the  next  two  “bullets”  and  the  stage  is  done.  Otherwise,  suppose  the  maximum  element  was  in  row  r ' 
Interchange  rows  r '  and  r+7.  This  is  the  pivoting  procedure.  To  make  the  permutation  a 
similarity  transformation,  also  interchange  columns  r '  and  r+7. 

For  /  =  r  +  2,  r  +  3,  .  .  .,7V,  compute  the  multiplier 


0,r 

*Lr+ 1  = - 


a, 


r+l,r 


Subtract  «,  ,  times  row  r+7  from  row  i.  To  make  the  elimination  a  similarity  transformation,  also  add 


«.  r+1  times  column  /'  to  column  r+7. 

A  total  of  N-2  such  stages  are  required. 

4.3.3  The  QR  Algorithm  for  Real  Hessenberg  Matrices 

The  basic  idea  behind  the  QR  algorithm  is  that  any  real  matrix  can  be  decomposed  in  the  form 

A  =  QR,  (17) 

where  Q  is  orthogonal  and  R  is  upper  triangular.  Now  consider  the  matrix  formed  by  writing  the  in  the 
opposite  order: 

A'  =  R  Q  .  (18) 

Since  Q  is  orthogonal,  equation  ( 1 7)  gives  R  —  Q1  •  A  .  Thus  equation  (18)  becomes 
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A'  =  Q'AQ.  (19) 

We  see  that  A'  is  an  orthogonal  transformation  of  A  . 

The  QR  algorithm  consists  of  a  sequence  of  orthogonal  transformations: 

A=QSRS,  (20) 

A.^=RSQS  (=Ql-As-Qs).  (21) 

The  QR  transformation  preserves  the  upper  Hessenberg  form  of  the  original  matrix  A  =  A]}  and  the 


workload  on  such  a  matrix  is  0(n 2  )  per  iteration  as  opposed  to  0(«3 )  on  a  general  matrix.  As 

s  00 ,  A  converges  to  a  form  where  the  eigenvalues  are  either  isolated  on  the  diagonal  or  are 

eigenvalues  of  a  2x2  sub-matrix  on  the  diagonal. 

In  order  to  accelerate  convergence,  we  deployed  the  technique  of  shifting:  If  k  is  any  constant,  then 

A  —  kl  has  eigen-values  \  —  k  .  If  we  decompose 

\ 

As  -  kj  =  QS  RS,  (22) 

so  that 


As+i  =  Rs-Q,+ksI 

=Ql-AQ, 


(23) 


then  we  verified  that  the  convergence  is  determined  by  the  ratio 

Aii 

here  A,  <  A( .  A  good  choice  for  the  shift  ks  would  maximize  the  rate  of  convergence. 


Any  real  matrix  can  be  triangularized  by  pre-multiplying  it  by  a  sequence  of  Householder  matrices 
(acting  on  the  first  column),  P2  (acting  on  the  second  column), . . . ,  Pn_t .  Thus  Q  =  Pn_x  •  •  •  P2  •  Pt ,  and 
the  first  row  of^  is  the  first  row  of  P]  since  is  an  (/ -l)x  (/  - 1)  identity  matrix  in  the  top  left-hand 


comer. 


The  Householder  matrix  is  determined  by  the  first  column  of(As  —  £J+)7)  •  (As  —  ksI ) ,  which 


has  the  form  [/?, ,  ,  rj ,  0, . . . ,  O]7  .where 

r  («»  -  a\  1  )(an-ln-l  ~an)~  an-\.n0n.n-\  ,  „  n 

P\  =«2l[ - +  a!2J 

°7\ 

<1\  =  a7\ \-a27  ~  a\  I  _  ( ann  ~  a\\  )  _  (^n-l.n-l  —  a\  I  )]  *  @4) 


Since  it  has  only  first  3  elements  nonzero,  the  matrix  P\  -  As-  TJ7  is  upper  Hessenberg  with  3  extra 


elements: 


x  x  x  x  x 

X  X  X  X  X 


-p?  = 


X  X 
X  X 


XXX 
XXX 
X  X 


X 


This  produces  a  matrix  with  the  3  extra 
xx  x  x  x  x  x 

XX  X  X  X  X  X 

X  X  X  X  X  X 

X  x  X  X  X  X  . 

X  X  X  X  X  X 

XXX 
X  X 


X  X 
X  X 
X  X 
X  X  . 

X  X 
X  X 
X  X 

elements  appearing  one  column  over: 


Proceeding  in  this  way  up  to  Pn_x ,  we  see  that  at  each  stage  the  Householder  matrix  Pr  has  a  vector  that  is 

nonzero  only  in  elements  r,  r+7,  and  r+2.  These  elements  are  determined  by  the  elements  r,  r+7,  and  r+2 
in  the  (r-1) st  column  of  the  current  matrix. 

In  summary,  to  carry  out  a  double  QR  step  we  construct  the  Householder  matrices, 

Pri  r  =  -1 .  For  Px  we  use/?,,  qx ,  and  rx  given  by  (24).  For  the  remaining  matrices,  pr , 
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qr ,  and  rr  are  determined  by  the  (r,  r-1),  (r+1,  r-1),  and  (r  +  2,  r-  1)  elements  of  the  current  matrix. 
Algorithm  description: 

There  are  two  possible  ways  of  terminating  the  iteration  for  an  eigenvalue.  First,  if  ann_x  becomes 

“negligible,”  then  an  n  is  an  eigenvalue.  We  can  then  delete  the  nth  row  and  column  of  the  matrix  and 

look  for  the  next  eigenvalue.  Alternatively,  anA  n_2  may  become  negligible.  In  this  case  the  eigenvalues 

of  the  2x2  matrix  in  the  lower  right-hand  comer  may  be  taken  to  be  eigenvalues.  We  delete  the  nth  and 
(n-1) st  rows  and  columns  of  the  matrix  and  continue. 

The  test  for  convergence  to  an  eigenvalue  is  combined  with  a  test  for  negligible  sub-diagonal 
elements  that  allows  splitting  of  the  matrix  into  sub-matrices.  We  find  the 
largest  /  such  that  at  M  is  negligible.  If  i  =  n,  we  have  found  a  single  eigenvalue.  If 

/  =  n- /,we  have  found  two  eigenvalues.  Otherwise  we  continue  the  iteration  on  the 
sub-matrix  in  rows  /  to  n. 

After  determining  /',  the  sub-matrix  in  rows  /  to  n  is  examined  to  see  if  the  product  of  any  two 

6 

consecutive  sub-diagonal  elements  is  small  enough  that  we  can  work  with  an  even  smaller  sub-matrix, 
starting  say  in  row  m.  We  start  with  m  =  n-2  and  decrement  it  down  to  i+1,  computing  p,  q,  and  r 
according  to  equations  (24)  with  1  replaced  by  m  and  2  by  m+1.  If  these  were  indeed  the  elements  of  the 
special  “first”  Householder  matrix  in  a  double  QR  step,  then  applying  the  Householder  matrix  would  lead 
to  nonzero  elements  in  positions  (m+l,m-l),  (m+2,m-l),  and  (m+2,m).  We  require  that  the  first  two  of 
these  elements  be  small  compared  with  the  local  diagonal  elements  am_{  m_, ,  amm  and  am+l  m+1 .  A 
satisfactory  approximate  criterion  is 

K.m-lltM  +  M)111  |/7|(|am+l,m+l|  +  |am,m|  +  |am-l,m-l|)-  <25) 

If  ten  iterations  occur  without  determining  an  eigenvalue,  the  usual  shifts  are  replaced  for  the  next  iteration 

by  shifts  defined  by 

K+K+i  =  1  -5  x  (|a„t(I_,  |  +  |tf„_i,n-2 1) 

•  (26) 

KK+\  —  |  |^w-l,n-2  |) 
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This  strategy  is  repeated  after  20  unsuccessful  iterations.  After  30  unsuccessful  iterations,  the  routine 
reports  failure. 

Because  our  system  carries  out  exact  arithmetic  operations,  we  are  able  to  set  the  tolerance  for  the 
condition  as  “negligible”,  which  means  using  an  <  Tolerance  and 


\am  „ 

m,n 

t-1 

kl 

K 

H 

D 

wd 

am+l,m+l 

1  + 

am 

m. 

m  | 

+  1 

a  _.=0 

<  Tolerance ,  instead  of  n'n  1  and 


+  H)  \p\  (|am+l,m+l|  +  |^m,m  |  +  |))  | P\  (|^m+l,m+l  |  +  |) 

when  deploying  floating  point  arithmetic.  The  results  would  be  more  precise  with  a  smaller  tolerance. 

4.3.4  The  Modules 

Our  program  calculates  the  eigenvalues  of  an  Nx  Aureal  general  Matrix  A.  The  algorithm  is  a  translated 
version  of  the  EISPACK  subprogram  RG.F. 

List  of  Routines: 

RG.F  calls  subroutines  BALANC,  ELMHES  and  HQR. 


Program  Overview  Flowchart: 
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BALANC 


void  balanc(vec_ZZ  &  vecTnum,  vecZZ  &  vecTden, matZZ  &  matPadicsl,int  &  nm,  int  &  n,int  & 
low,int  &  igh,vec_ZZ  &  scaleNum,  vecZZ  &  scaleDen); 

Function  Description:  It  balances  a  real  matrix  and  isolates  eigenvalues  whenever  possible. 

INPUT: 

nm  must  be  set  to  the  row  dimension  of  two-dimensional  array  parameters  as  declared  in  the  calling 
program  dimension  statement. 
n  is  the  order  of  the  matrix. 

vecTnum  vecTden  and  matPadicsl  contain  the  input  matrix  to  be  balanced. 

OUTPUT: 

vecTnum  vecTden  and  matPadicsl  contain  the  balanced  matrix. 

low  and  igh  are  two  integers  such  that  A[i,j]  is  equal  to  zero  if  ( 1 )  i  is  greater  than  j  and  (2)j=l . low-1  or 

i=igh+ 1 

scale  contains  information  determining  the  permutations  and  scaling  factors  used. 

ELMHES 

void  elmhes(vec_ZZ  &  vecTnum,  vec  ZZ  &  vecTden,  matZZ  &  matPadicsl, int  &  nm,  int  &  n,int  & 
low, int  &  igh  ,int  &  inte  ); 

Function  Description:  Given  a  real  general  matrix,  it  reduces  a  sub-matrix  situated  in  rows  and  columns 
low  through  igh  to  upper  Hessenberg  form  by  stabilized  elementary  similarity  transformations. 

INPUT: 

nm  must  be  set  to  the  row  dimension  of  two-dimensional  array  parameters  as  declared  in  the  calling 
program  dimension  statement. 
n  is  the  order  of  the  matrix. 

low  and  igh  are  integers  determined  by  Balanc.  If  Balanc  has  not  been  used,  set  /ow=7,  igh=n. 
vecTnum  vecTden  and  matPadicsl  contain  the  input  matrix. 

OUTPUT: 

vecTnum  vecTden  and  matPadicsl  contain  the  Hessenberg  matrix.  The  multipliers,  which  were  used  in  the 
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reduction,  are  stored  in  the  remaining  triangle  under  the  Hessenberg  matrix. 

inte  contains  information  on  the  rows  and  columns  interchanged  in  the  reduction.  Only  elements  low 
through  igh  are  used. 

HQR 

void  hqr(int  &  nm,  int  &  n,  int  &  low, int  &  igh,vec_ZZ  &  vecTnum,  vec  ZZ  &  vecTden, matZZ  & 
matPadics  1  ,vec_ZZ  &  wrNum.vecZZ  &  wrDen,vec_ZZ  &  wiNum,vec_ZZ  &  wiDen,int  &  ierr); 

Function  Description:  It  finds  the  eigenvalues  of  a  real  upper  Hessenberg  matrix  by  the  QR  method. 

INPUT: 

nm  must  be  set  to  the  row  dimension  of  two-dimensional  array  parameters  as  declared  in  the  calling 
program  dimension  statement. 
n  is  the  order  of  the  matrix. 

low  and  igh  are  integers  determined  by  Balanc.  If  Balanc  has  not  been  used,  set  low=l,  igh=n. 
vecTnum  vecTden  and  matPadics  1  contain  the  upper  Hessenberg  matrix.  Information  about  the 
transformations  used  in  the  reduction  to  Hessenberg  form  by  Elmhes,  if  performed,  is  stored  in  the 
remaining  triangle  under  the  Hessenberg  matrix. 

OUTPUT: 

vecTnum  vecTden  and  matPadics  1  have  been  destroyed.  Therefore,  it  must  be  saved  before  calling  hqr,  if 
subsequent  calculation  and  back  transformation  of  eigenvectors  is  to  be  performed. 
wrNum  wrDen  and  wiNum  wrDen  contain  the  real  and  imaginary  parts,  respectively,  of  the  eigenvalues. 
The  eigenvalues  are  unordered  except  that  complex  conjugate  pairs  of  values  appear  consecutively  with  the 
eigenvalue  having  the  positive  imaginary  part  first.  If  an  error  exit  is  made,  the  eigenvalues  should  be 
correct  for  indices  ierr+ 

ierr  is  set  to  zero  for  normal  return,  or  set  to  j  if  the  limit  of  30*n  iterations  is  exhausted  while  they'-th 
eigenvalue  is  being  sought. 

The  matrix  and  related  parameters  are  stored  as  p-adic  sequences  in  matZZ  and  fractional  number. 
For  the  rational  data  structure,  the  numerator  and  denominator  of  the  rational  are  stored  separately  in  array 
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vecTnum  and  vecTden  as  arbitrary  length  integers.  This  maintains  the  precision  of  the  numbers  during  the 
calculation.  During  computational  process  all  rational  numbers  will  keep  their  fractional  data  type.  The 
floating-point  data  type  was  never  used,  to  keep  the  calculation  free  of  round-off/truncation  errors. 


Note  the  Error  Code  output:  if  normal  return,  ierr  =  -/;  if  Error  Code  >  0,  it  indicates  that  more  than  30 
iterations  of  a  subroutine  were  required  to  determine  an  eigenvalue.  In  this  case,  the  subroutine  terminated 
after  30  iterations. 

The  eigenvectors  are  outputted  as  a  square  NXN  matrix  whose  entries  correspond  to  the  eigenvalues  as 
follows:  If  the  i-th  eigenvalue  is  real,  the  i-th  COLUMN  of  the  eigenvector  matrix  contains  the 
corresponding  eigenvector;  If  the  i-th  eigenvalue  is  complex  with  positive  imaginary  part,  COLUMNS  i 
and  (i+1)  of  the  eigenvector  matrix  contain  the  real  and  imaginary  parts  of  the  corresponding  eigenvector. 
4.3.5  Result  Analyses 

We  use  the  following  example  to  show  the  advantages  of  our  programs  in  terms  of  exactness. 


The  input  matrix 
Result: 


4 

111 

5 

111 

1 

111 


_14 

111 

n_ 

in 

o 


—  0 


0 

_2_ 

111 


-  has  eigenvalues 


4=777  4=777  4=777 


111 


111,  and 


THE  BALANCED  MATRIX  is: 

2/111  0/1  -1/111 
0/1  13/111-5/111 
0/1  14/111  -4/111 

THE  UPPER  HESSENBERG  MATRIX  is: 
2/111  0/1  -1/111 
0/1  13/111  -5/111 
0/1  14/111  -4/111 


Eigenvalue  1: 
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2/111  +  jO/l 
Eigenvalue  2: 

4403202284312065429037843552686838343371 5773553 1 7037774709535840028 1 485847264778 
3885308324 1 1 739293057334922 1 73 1 3050558 1 325899767 1 2796782 1 74639787 1195589181 88959 
22170272007433966278435458954982246587849564662126277527673707724488387314845240 
7316611 94486786 1 1 7403443067224453 1 9604 1 004464557 1 1 872024503408794989428024999557 
72367374331500429977063062653905440955495027034700110517834687322583943548114289 
20487384484543724036357178869514787915773160244610311077901974476456747955678844 
300477672788033910258324348335702228394885570985506 1 30783 1 859687008 1 849669675244 
73 1 356390848297884962875487623802 186511 698999 1081448516161 7327888699986994382 1 6 1 
0883 1 02477 1 732 1 665797227094 1 295286467695635871436909 158550227704979 1 4672443 15516 
3328839 13801 46088376 1 5600860829960703893290772819960569052 1018011 7357678 1 5276798 
5909433235933007948 1 4674596628233 1 8370479626325472496352345432847465969088296 1 04 
28373211 047627295025895 1 1 493920092 1 720 177457521 6626966674525846 1 34285667095275 1 2 
18032601247038636424675563928171582647205470537422467^488808855048852276404241 
1 56655 1 1 82469544603 1 0972 1 1 30 1380086764455 1 57 1 864677 1 295527323482599 1 2807932883 1 3 
3873 1 00486 1 252085 1 23399344975 1 8382644559632777422759647923 1 00 1 87946764797 1 989059 
857053243372044294 1 0057744644 1331511845031 8432 1 6 1 37246487876732 135061 74227469332 
53386680227576456 1 34520568 1 478160083845750294554867 1 747058388995544208967530585 1 
35956 1 788 1 9028560493 1426578 1 4827072 1 0307 1 06 1 0036384625 1 24036303 1 5276086657439555 
78853367947057069029770 1 6633360287989 1 1 90733 1 34 1 52298 1 50652563067540320404704 1 45 
78228990 1 480576650952 1 08 1 73473689 1 569672971454632686222373923381428445230260877 1 
59737283720351212204287788577980070266208143665421741015899903767280098068319985 
1 600 1 2923046320655 1 473746264382582 1 745396790 1 3077622765637238847662928796 1 434674 
14325788357 161342537827867288662463808357 1 7826393230159659806355454822766 1 025957 
24439 1 99876 1 0032 1 325798 1 96 1 26737664 1 246 1 640 1 696776542557373577409745580597249478 
970873496585370487704703 1 03920 1 1 65986366 1 90952 1 7844707145734563280289345 1 3404 1 99 
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30608111728960423077720099236589348239679428609695690342515610380878344966858032 
041402236574245 1 0 1 1429838358 1 1 70479103307987222089797995959058963603843466968 1 75 
732226750671036838008333496440823855 1 3285097795 1 72673838 1 886443521/8145920243782 
43286698659646628692 1 929849825 1 35680686677328242730575092595445698763760 103021 29 
6646704 1 03889430562 1 45990303922254 1 390787332954482 1 005063098098383936359732 13081 
660799689 1 84848285 14581 68564947780 1 3 1 2200555078303900637798083 1 852 1 929 1 9063053 1 1 
236849018099209545429677 1 5335547269399378597 1 587063245 158897263667 1 6 1 83278863345 
256 1 79840992 1 87 1 7269 1 59098033 1 069 1 98063462 1 564656396848965373 1 8950 1375452 1 469 1 32 
99324 1 38603 1 98 1 7695053 1 00379 1 82097 1 646699862 1 23473 1 23343263389 1 9745381 8707 1 26769 
745386648508485590488554440598748219822363174558644281 137758767133886462251471 1 1 
224434609560387090 1 7 1 49353546276946 114981 86659 1 2 1 454975274725335 14281 4554245334 1 
04 1 1 2587479588894604427 1 056645232990004396 1 26043 1 85057473453656 1 973069376639 1 724 
30457710640595 101071561 0420060945983 1 96065659384 176111 33067664 1 88736335746269603 
472 1 70 1 8237629 1 4464776940733 1 04762872965750 1 1 2072 1717001 97982750474465 1812511147 
77306329148716635243562438340770545006545492735143929314507702854422264288277652 
90208206058948 1 96325 1 08270367 1 0726 1 2652587464 1 33527557322724373596 1 832573747676 1 
2982 1 2 1 74058400 1 647370732862200667749943 1 97634089330660078587 1 48064376795 1 494 1 67 
23211 2089 1 963350047736 1 3 1288332789926426848597953 14715321 0920589837763 1 543296859 
69 1 2480847598459759 1 3065382888482658 1 209054565890659 1 3966 1 059767668782626382 1 894 
53452 1 56 1 495740762928 1 656684 1 70683660968340705244222043222 1 45344 1 998099696529593 
39804726317722842121189777843593316564989925883254903846988982682944839850835487 
592832 1 52745908860488 1 83994094908260 1 0 1 45312968402 1 64085993658489747629840 1 89356 
2082 1 96536503 1 6384656844 1 346 1 095636229 1 6767 11621 36793 1 03060 1 98 1 8722940292786920 1 
3184121 2474 1 83 1 9722 166361 30499323067074988797 17831 08960972536868728833381 6007879 
729665 1 651497438985974538778305 1481 788526044206996095620 191371 090402633048892058 
02607758563696 1 1 8773 104861 8300786 1 868344 1 93 1 6592 1 7660545 1 158923 1 682963587 1 775756 
3963839232495403650640979 17511 5354636 1 9574 1 5 1 37723745351 90926748959222623342543 1 
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63 1 78856539022 1 375346032328444479 1 9243 1 4 1 7528734543883494 13852459 13831611 9468484 
4268 1 807604340724694 1 07445 1 2524722 1 1 004396524945825223836 1 877844 1 773705 1 74597789 
537202756767478525580308 1 25564 1 579869865903673282257062489 1 802079542434 148013549 
68779901 1 16559585476875 1548809514622895473538432573440+ jO/1 
Eigenvalue  3: 

220 1 59920487203306097293653549728073663236506982595 1381 7893848302 144612882181271 
8 1 96368 1 686863932079703089863077985934 192 1 153761933520222468486 1 977 1 4023 1 5365203 
02874463 5  8499 1249811 43806703 12861891 225404 1900115927171 863702899 1 8572947 10541184 
64586069381 502 1 39974761 1 383304320705877545237 1 7307653392222573546440136938850449 
5897 1 704 1 97 1 87560775968790984 1 42809736084680874 1 1 053479795 1 633588352405735278906 
6088 1 0965965592342624702570121 53840396037365030506255562 1 3 1 703496498864982608338 
4 1 1 5298 1 8 1330443 1 9571462565370 1 2 1 254533750575573248365354677 1 1 4 13334957336942985 
8794359643309566293201 797697278388 1 96754239945467439424237642522802 1 49020376 1 529 
366350939 1 837 1 65776425857 1 9 1 07725830 1 54389738639334658334392992695522028768 1 7456 
04 1 1 37772455192625922 1 740929 1 32688 1 775659565835 1 248953405 1 292 1219631 0853548498 1 3 
6826589755704 1 6 1 660289484674 1 90759 1163516146381 478930742900 1 3224 1 84073570764 1 248 
5668656022977608608 1 33026996873688570082 1 270693  52442895 1 5  8445964308 1711 52876029 1 
346 1 3226797874339342856966503 147183019135671 057200674 1 729 1 647740285392099 1 796339 
0234549 1 1 030146237008403 1 85072764 1592101 576447729 1 26556268566844668 1363752 1 08057 
3 1 48682 1 54063705 1 90 1 344460396436 178891 03687403 1 354250 1 454633602387 1 99082860859 1 8 
75276729547 13834782 125411 220385 1 1 8146589926483563726049 1 354823467692500554304646 
3942 1 74650 1 969855759 1 936439121 5270664 1 3 1 363882066740571 805548 1 7474658 127411261 89 
0934578 1 85959970 1 4380773 1 257990209390472072338 1512641531 97340892326590534 1587131 
75787656885711407787753940590745006776467952245982817538748248391777244502585410 
4906978864575308205752802595 1 5723859360 1 703234755607948369597799395069227778722 1 
440 1 1 22209 1 75782090699984483245659336658403653794 1 398 1 337344 1 9465014920385202738 
88 1 964258 1 9 1 836964895030407235979832396253242900467885828448984734595373686 14312 
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0924 1 27386 1 742076335285465328 151003361 827240627227895455333385 1 6707209295 1 972598 
0034477733702342 1 445523557 1 5 1 08092 1 853535927872 1 9347820265 1 863704794563347086 1 88 
574007903 1436068649023809854 1 778 1 00536 1330741981 68310502 1 48244979572992273947 1 43 
222676067767 1 844772 171433363791 69470 1 22694567993 1 0393900 1 42 1 40943332 1 894 1 3289089 
68801516688438755319043299165764361209017387501371042592592567782859133025205788 
5282300 1 3904272779343 1 557385 1 930785667085574736853 1 246603 1 69 1 694076/814592263309 
9365773026644929997 1 593330823809 1 8674 1 8660085888 1 87466 1 2819672246699066309636062 
2 1 7202 1 1 005732 1 334423445 1 564874377580683592 1 53202944772296279 158328048 1 528947805 
78345767820098833170876979938995887012656077920947730207621842805625983820178092 
78206427854577300096023469 1 649 1 74566 161530910851 320 1 54 1 8900696804 1 777206090096 1 3 
4650 1 99390396829452580582452080428530937652342236548680776348768266 1 567590826543 
755066459005964722394358266598807 1038331 3290956 1 30593800720766371 4 1 532627850 1 865 
5074558526803660762490 1 56959469337227046996963 1 877 1 029 1 5404 1 68408057 1 97654599726 
05962969317420674602728676940439978848677804058882538447095453898255933846073099 
8686940699 1391004161 856875648 1 39 1 049097092444 1 1 7002654984224558 15012357251 799988 
5937640595879335925542625489577 1 362234346579 1 1 586669747202498 1 524795889376939972 
70276905 15838539151 9794 1 67525627 1 2642365725 1 14825932936730488 11088655410151 43074 
5738973 1 909526247922765606930608670845855 120581 7746780782893556897 1 54 1 0733404948 
24494759632233888 1 2346 1 0692 1 2756586989549450779 1 1 1 9393335 1 69243030924 1 69 1 6680324 
57 1 379432141 381451 9084460 1 690 1 04744079 138659406926389 11435311 54092536 1 06322 1 3726 
58742935330058262473141 80712725694536320283 13855933450703736171437915465908320 1 7 
3568 1 637503 1 3593889 1 4500763300885276237 1 798 1 9806735306477 1 78 1 639802662564270340 1 
880043955638 1 3607 1 2748 1 96955 1 9509690546437286302 1 5539358008397727 1 74933 1 0 1 446074 
9529433644 1 059088 1 4289762552262 1 44670 1 292985803359846 1 103357490652247 1 88664 1 763 1 
8371 795 1 65586826574 11494151 980934 1 1 94945675974308 1091533195395561191 7235353982 1 7 
68841 134298608521 1 18123799908426971871983775197035005579864463823225649200964063 
951191 952584 1 26494467445400 1 3036760886432482837606425265297870582675060966984896 
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605113820191825102410519529017591044331 5374034 163105288331 7292036 1 04404338206446 
8 1 1 9452733 1 480535 1 6625902 1 74967 10177121 924494902 1 925 1 97033733943 1 485726709989654 
96942207328995862 1 39063569994252025276440535290 1 7987463 1 6708336465053467 1 89 1 8227 
0265 1 289483605420634326577978497485950488 1 632914746996460938639770267623948 1 8930 
83776888573483781601940843849175645606233866292110823670307664350086514163466508 
9935503789827568920 1 553037 1 74 1 22448467622 1 3698386835242526 1 9204 17856760 137613786 
8384021496968048835052292594 1 6966 1 1 133249854543 1 2552459 1  +  jO/1 . 

As  far  as  exactness  is  concerned,  as  we  can  see  from  the  above  example,  our  system  can  obtain  much  more 
accurate  results,  since  floating-point  arithmetic  operations  can  only  provide  the  precision  limited  by  the 
precision  of  the  data  type.  Our  system  performs  much  more  accurate  results  and  can  be  trusted  in  the 
process  of  exact  computation. 

5.  Conclusions  and  Discussion 

✓ 

Despite  the  advantages  of  finite  segment  p-adic  arithmetic,  there  is  a  problem  need  to  be  solved — the 
estimation  of  an  efficient  length  M  of  the  p-adic  expansion  in  arithmetic  operations.  Currently,  it  cannot 
give  the  sufficient  number  of  digits  for  p-adic  sequence  for  exact  computation  in  some  situation.  As  we 

concluded  that  the  maximum  length  of  the  p-adic  expansion’s  periodical  part  of  (tf  ±  b )  is  LCM ( n ,  m) ; 
in  multiplication,  the  maximum  length  of  periodic  part  of  the  product  of  ( a  x  b)  is 

LCM{m,ri)x(pGCD{m'n)  - 1) . 

In  order  to  compute  matrix  exponential  with  JCF,  we  performed  eigenvalue  analysis.  The  routines  we 
used  are  translations  of  EISPACK:  The  matrix  is  first  balanced;  Orthogonal  similarity  transformations  are 
used  to  reduce  the  balanced  matrix  to  a  real  upper  Hessenberg  matrix;  The  implicit  double-shifted  QR 
algorithm  is  used  to  compute  the  eigenvalues  and  eigenvectors  of  this  Hessenberg  matrix. 

Thus,  to  determine  the  efficient  length  M  of  the  p-adic  expansion  for  a  rational  number  in  the  matrix 
operation  becomes  an  important  problem.  In  particular,  the  algorithm  in  eigensystems  may  use  iterative 
arithmetic  operations.  This  made  the  situation  even  more  complicated.  From  our  experience  of 
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implementing  the  eigensystem  in  p- adic  field,  the  mechanism  for  setting  the  length  M  is  our  biggest 
concern,  because  running  time  is  related  to  it,  a  large  M  will  increase  the  runtime  tremendously.  When 
dealing  with  iterative  arithmetic  operations,  the  required  length  of  the  p-adic  sequence  will  keep  on 
increasing  during  the  computation. 

One  possible  solution  is  to  create  a  gradual  increasing  “AT.  As  we  know,  the  corresponding  fraction 
numbers  will  be  completely  different  for  the  same  p-adic  sequence  with  different  lengths.  If  we  can  predict 
the  range  of  the  approximate  eigenvalue  based  on  previous  results,  we  would  be  able  to  adjust  the  “Af” 
when  an  abnormal  result  occurs.  To  be  precise,  we  can  set  a  judgmental  mechanism  after  the  operation,  if 
an  abnormal  result  occurs,  redo  the  operation  with  an  larger  “A f\  at  the  same  time,  reduce  the  p-adic 
expansion  of  the  products  to  its  minimum  size. 
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